/*
    Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    This is the old code that was used to calculate the interpolation
    file used when projecting particles. It's based on the very first
    iteration of grids, so is totally outdated. If we ever need to
    resurrect it, it'll require a serious amount of figuring out what
    we were trying to do. */

#include <iostream>
#include <cmath>
#include "grid.h"
#include "optical.h"
#include "sphgrid.h"
#include "sphparticle.h"
#include "misc.h"
#include "counter.h"
#include "integrator2.h"

using namespace mcrx;
using namespace std;

typedef particle<nbody_data> sphparticle;
typedef grid_cell<nbody_data_cell> T_cell;
typedef nbody_data_grid sphgrid;
typedef nbody_data_grid::T_cell_tracker T_cell_tracker;

class sphintegrator : public integrator2 {
public:
  class yintegrator : public integrator2 {
  public:
    class zintegrator : public integrator2 {
    public:
      sphparticle* p;
      double x,y;
      double integrand(double); // The integrand
    };
    double integrand(double);
    zintegrator zi;
    double zlo,zhi;
  };
  double integrand(double);
  yintegrator yi;
  double ylo,yhi;
  double project(sphparticle*,T_cell*);
};

double sphintegrator::project(sphparticle* p, 
			      T_cell* c)
{
  // Avoid integrating where the kernel is identically zero
  // so check what the overlapping cube is
  double pminx=p->position()[0]-2*p->radius();
  double pminy=p->position()[1]-2*p->radius();
  double pminz=p->position()[2]-2*p->radius();

  double pmaxx=p->position()[0]+2*p->radius();
  double pmaxy=p->position()[1]+2*p->radius();
  double pmaxz=p->position()[2]+2*p->radius();

  vec3d cellmin=c->getmin();
  vec3d cellmax=c->getmax();
  int notwithin=0;
    
  if(cellmin[0]>pminx) {
    pminx=cellmin[0];
    notwithin=1;
  }
  if(cellmin[1]>pminy) {
    pminy=cellmin[1];
    notwithin=1;
  }
  if(cellmin[2]>pminz) {
    pminz=cellmin[2];
    notwithin=1;
  }
  
  if(cellmax[0]<pmaxx) {
    pmaxx=cellmax[0];
    notwithin=1;
  }
  if(cellmax[1]<pmaxy) {
    pmaxy=cellmax[1];
    notwithin=1;
  }
  if(cellmax[2]<pmaxz) {
    pmaxz=cellmax[2];
    notwithin=1;
  }
  
  if( notwithin && (pmaxx>pminx) && (pmaxy>pminy) && (pmaxz>pminz) ) {

    // We have a nonzero region where kernel and cell overlap - integrate
    ylo=pminy;
    yhi=pmaxy;
    yi.zlo=pminz;
    yi.zhi=pmaxz;
    yi.zi.p=p;

    // Because of how the integrator works we don't want the identically
    // zero region in the beginning since then the step size is increased and
    // overshoots when the kernel goes nonzero leading to a huge penalty
    // so, we divide the integral in two 
    if(pminx>= p->position()[0])
      return integrate(pminx,pmaxx);
    else if(pmaxx<= p->position()[0])
      return -integrate(pmaxx,pminx);
    else
      return -integrate((double)p->position()[0],pminx)+
	integrate((double)p->position()[0],pmaxx);
    
  }
  else if( notwithin )
    return 0;
  else
    return 1;
    
}

double sphintegrator::integrand(double x)
{
  yi.zi.x=x;

    if(ylo>=yi.zi.p->position()[1])
      return yi.integrate(ylo,yhi);
    else if(yhi<=yi.zi.p->position()[1])
      return -yi.integrate(yhi,ylo);
    else
      return -yi.integrate((double)yi.zi.p->position()[1],ylo)+
	yi.integrate((double)yi.zi.p->position()[1],yhi);
  
}

double sphintegrator::yintegrator::integrand(double y)
{
  zi.y=y;


    if(zlo>=zi.p->position()[2])
      return zi.integrate(zlo,zhi);
    else if(zhi<=zi.p->position()[2])
      return -zi.integrate(zhi,zlo);
    else
      return -zi.integrate((double)zi.p->position()[2],zlo)+
	zi.integrate((double)zi.p->position()[2],zhi);
}

double sphintegrator::yintegrator::zintegrator::integrand(double z)
{
  double px=p->position()[0];
  double py=p->position()[1];
  double pz=p->position()[2];

  // Here's a function of the coordinates x,y,z and the sphparticle p
  double r=sqrt( (x-px)*(x-px) + (y-py)*(y-py) + (z-pz)*(z-pz) );

  return p->kernel(r);
}

void makeinterpolationdata();
void testinterpolationdata();
sphintegrator sphint;

int main(int, char**)
{
  sphint.setaccuracy(1e-4,1e-4);
  sphint.yi.setaccuracy(1e-4,1e-4);
  sphint.yi.zi.setaccuracy(1e-4,1e-4);

  //makeinterpolationdata();
  testinterpolationdata();
}

void testinterpolationdata()
{
  ifstream inf("/home/sunrise/dust/SNAP010/SNAP010-550.sph");
  double q;
  inf >> q >> q >> q;
  sphparticle p(0,vec3d(0.5,0,0), vec3d(0,0,0), 0.25, 
		nbody_data());

  sphgrid g(vec3d(-1.5,-1.5,-1.5), vec3d(1.5,1.5,1.5), vec3i(1,1,1));
  T_cell_tracker c =g.locate(vec3d(0,0,0),0,false);
  c.cell()->refine(vec3i(2,2,2));
  //c =g.locate(vec3<double>(0,0,0));
  //->refine(c,vec3<int>(2,2,2));

  while(inf.good()) {
    inf >> p;
    g.init();
    while(c=g.enumerate()) {
      double a1=p.project2(c);
      double a2=sphint.project(&p,c);
      double r=mag(c->getcenter()-p.position());
      double a3=p.kernel(r)*c->volume();
      if(a2!=0)
	cout << a1 << '\t' << a2 << '\t' << a3 << '\t' << p.radius() << '\t' 
	     << r << '\t' << c->getsize()[0] << '\n';
    }
    }
}

void makeinterpolationdata()
{
  int nx=13;
  int ny=13;
  int nz=13;
  int nh=20;

  cout << "4\n" << nx << ' ' << ny << ' ' << nz << ' ' << nh << "\n";

  double xx[nx];
  double yy[ny];
  double zz[nz];
  double hh[nh];

  int i=1;
  xx[0]=-1e30;
  xx[nx-1]=1e30;
  for(double x=-2;x<=2.01;x+=4./(nx-3)) {
    xx[i++]=x;
  }

  i=1;
  yy[0]=-1e30;
  yy[ny-1]=1e30;
  for(double y=-2;y<=2.01;y+=4./(ny-3)) {
    yy[i++]=y;
  }

  i=1;
  zz[0]=-1e30;
  zz[nz-1]=1e30;
  for(double z=-2;z<=2.01;z+=4./(nz-3)) {
    zz[i++]=z;
  }

  i=1;
  hh[0]=0.0001;
  for(int l=1;l<nh;l++) {    
    hh[l]=0.01*pow(10.,4.*(l-1)/(nh-2));
  }

  for(int i=0;i<nx;i++)
    cout << xx[i] << ' ';
  cout << '\n';
  for(int j=0;j<ny;j++)
    cout << yy[j] << ' ';
  cout << '\n';
  for(int k=0;k<nz;k++)
    cout << zz[k] << ' ';
  cout << '\n';
  for(int l=0;l<nh;l++)
    cout << hh[l] << ' ';
  cout << "\n\n";
    
  sphgrid g(-0.5,0.5,1,-0.5,0.5,1,-0.5,0.5,1);


  sphparticle p(0,vec3d(0,0.5,0), vec3d(0,0,0), 0.001,
		nbody_data(1,vec3d(0,0,0),0,0,0,0,vec3d(0,0,0),0,0,0,0,0));
  T_cell* c =g.locate(vec3d(0.));
  //cout << sphint.project(&p,c) << '\n';

  double y=0;
  counter cc(10);
  for(int l=7;l<nh;l++)
    for(int k=1;k<nz;k++)
      for(int j=1;j<ny;j++) {
	for(int i=6;i<nx;i++) {
	  //cerr << l << ' ' << k << ' ' << j << "    \r";
	  //cc.incr();
	  p.h=hh[l];
	  p.pos[0]=xx[i]*p.h+0.5;
	  p.pos[1]=yy[j]*p.h+0.5;
	  p.pos[2]=zz[k]*p.h+0.5;
	  if(p.position()[0]<0)
	    p.position()[0]=0;
	  if(p.position()[1]<0)
	    p.position()[1]=0;
	  if(p.position()[2]<0)
	    p.position()[2]=0;
	  double a;
	  a = sphint.project(&p,c);
	  cout << a << ' ';
	  cerr << xx[i] << ' ' << yy[j] << ' ' << zz[k] << ' ' <<p.position() << '\t' << p.h << '\n';
	}
	cout << '\n';
      }
  /*
 
  grid<double> g(-5,5,10,-5,5,10,-5,5,10);
  grid_cell<double>* c;
  c=g.locate(matrix<double>(0.5,0.5,0.5));
  c->refine(matrix<int>(10.,10.,10.));
  c=g.locate(matrix<double>(-.5,0.5,0.5));
  c->refine(matrix<int>(10.,10.,10.));
  
  c=g.locate(matrix<double>(0.05,0.55,0.55));
  c->refine(matrix<int>(10.,10.,10.));
  c=g.locate(matrix<double>(-0.05,0.55,0.55));
  c->refine(matrix<int>(10.,10.,10.));
  
  matrix<double> p(0.005,0.555,0.555);

  grid_cell<double>* d;
  grid_cell<double>* dd;
  int t1,t2,t3,t4;
  for(int i=0;i<10;i++) {
    c=g.locate(p);
    t1=realcc();
    d=c->neighbour(p,0,1);
    t2=realcc();
    dd=c->neighbour(p,0,-1);
    t3=realcc();
    t4=realcc();
    }
  cout << t2-t1 << '\t' << t3-t2 << '\n';
  cout << c->getmin() << c->getmax() << '\n';
  cout << d->getmin() << d->getmax() << '\n';
  cout << dd->getmin() << dd->getmax() << '\n';
  */  
  /*  
  const int n=300;
  cout << n+1 << '\n';
  for(double i=0;i<=1;i+=1./n)
    cout << i << '\t' << cos(2*pi*i) << '\n';
  */

  //linear_interpolator i("sin.inp");
  //cout << i.interpolate(0.654);
}
